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Abstract 

We study a two-stage mixed-integer linear program (MILP) with more than 1 million binary 
variables in the second stage. We develop a two-level approach by constructing a semi-coarse 
model (coarsened with respect to variables) and a coarse model (coarsened with respect to 
both variables and constraints). We coarsen binary variables by selecting a small number of 
pre-specified daily on/off profiles. We aggregate constraints by partitioning them into groups 
and summing over each group. With an appropriate choice of coarsened profiles, the semi- 
coarse model is guaranteed to find a feasible solution of the original problem and hence provides 
an upper bound on the optimal solution. We show that solving a sequence of coarse models 
converges to the same upper bound with proven finite steps. This is achieved by adding violated 
constraints to coarse models until all constraints in the semi-coarse model are satisfied. We 
demonstrate the effectiveness of our approach in cogeneration for buildings. The coarsened 
models allow us to obtain good approximate solutions at a fraction of the time required by solving 
the original problem. Extensive numerical experiments show that the two-level approach scales 
to large problems that are beyond the capacity of state-of-the-art commercial MILP solvers. 

Keywords: Coarsened models, distributed generation, large-scale problems, two-level ap¬ 
proach, multi-period planning, resource and cost allocation, two-stage mixed-integer programs. 
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1 Problem Definition and Motivation 

We consider a hierarchical two-stage mixed-integer linear program (MILP) with integer variables 
in both the first and second stages. We are particularly interested in applications where the first- 
stage integer variables model design or purchasing decisions and the second-stage variables model 
operational decisions over a long time horizon (e.g., hourly operations over a decadal horizon). 
The goal is to take operational constraints into account when making a capital investment decision. 
Thus, our model is complicated by the fact that second-stage variables include both binary (on/off) 
decisions and continuous variables that model operational settings. Examples include the design 
of cogeneration units for commercial buildings subject to operational conditions (Siddiqui et ah, 
2005; Stadler et ah, 2009; Ren and Gao, 2010; Pruitt et ah, 2013, 2014) and transmission network 
expansion subject to unit commitment constraints (Alguacil et ah, 2003; Hedman et ah, 2010; 
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Munoz et al., 2013, 2014). Models of this class can involve hundreds of thousands or even millions 
of binary variables and are beyond the scope of today’s state-of-the-art solvers. 

We consider a two-stage MILP with m first-stage variables and three sets of N second-stage 
variables. The first-stage variables are y G {0,1}”^. We have three classes of second-stage variables: 
(1) on/off decisions, x G {0,1}'^; (2) operational settings v G that are switched on or off by x, 
that is. Lx < v <Ux for finite bounds L <U] and (3) other second-stage variables 0 < u; G . 
The MILP model is described by 


minimize a^y + x + c^v + d'^w 

y,x,v,w 

subject to Ay + Bx + Cv -|- Dw < f 

yG{0,ir 

X G {0,1}^ 

V G Lx < V < Ux 
w G w > 0, 


( 1 . 1 ) 


where a G b,c,d G M'^, A G and B,C,D G We assume that m N and 

m M; that is, the number of first-stage variables is much smaller than the number of second- 
stage binary variables and coupling constraints. 

In addition to a large number of binary and continuous variables in the second stage, the chal¬ 
lenges of model (1.1) also arise from coupling constraints. We do not assume any sparsity structure 
in matrices B, C, and D. This is in contrast to the arrow-type sparsity structure that arises, for ex¬ 
ample, in stochastic programs (Birge and Louveaux, 2011). Therefore, fixing the first-stage variable 
y does not decompose (1.1) into scenario-based subproblems. Moreover, the coupling constraints 
in (1.1) are not amenable to decomposition methods such as Lagrangian relaxation (Geoffrion, 
1974; Fisher, 1981). The reason is that the number of coupling constraints, M, is of the same 
order as the number of variables, N. Dualizing all coupling constraints results in a large num¬ 
ber of dual variables; hence, Lagrangian relaxation of (1.1) does not yield efficient decomposition 
methods (Fisher, 1981). 

Problem (1.1) arises in several applications. In particular, our work is motivated by the co¬ 
generation problem with renewable energy for commercial buildings. In this case, the first-stage 
design involves investment decisions for cogeneration units such as fuel cells, solar panels, and bat¬ 
tery storage. The second-stage problem aims at optimal on/off hourly operation that takes into 
account technology specifics such as minimum/maximum power generation. Our goal is to include 
operational constraints in the design of cogeneration. 

One of our objectives is to find the optimal first-stage solution of (1.1). However, the two-stage 
MILP with a large number of binary variables at the second stage is beyond the scope of state-of- 
the-art commercial MIP solvers. For example, a typical cogeneration model with a ten-year horizon 
results in 1.05 million binary variables (3650 days x 24 hours x 12 units). On the other hand, 
a naive approach that solves (1.1) with a short horizon at the second stage provides first-stage 
designs that are suboptimal for a long horizon problem. The reason is that short horizon problems 
do not take into account coupling constraints over a long horizon. Moreover, the problem data 
of short-horizon problems are not representative of long-horizon problems, resulting in suboptimal 
solutions. 

We develop a two-level approach that coarsens the hourly on/off variables to daily operation 
profiles. Since the profile representation yields a model with many fewer variables, we refer to this 
step as the primal (variable) coarsening. The resulting semi-coarse model still contains the same 
order of constraints as in (1.1). We reduce the number of constraints by partitioning them into 
groups that are of the same size as the profiles and summing over each group. This aggregation of 
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constraints results in a relaxed problem whose solutions may not be feasible for the original MILP 
model. We include the violated constraints, re-solve the coarsened MILP model, and repeat this 
process until all constraints are satisfied. We refer to the aggregation of constraints as the dual 
(constraint) coarsening and the resulting MILP as the coarse model. 

1.1 Literature Review 

The idea of aggregating variables and constraints to build approximate optimization models is not 
new. Zipkin studies the effect of variable aggregation and row aggregation for linear programs 
in (Zipkin, 1980b, a). This framework is extended to stochastic linear programs by Birge (1985) 
and by Clay and Grossmann (1997). For integer programs, a large amount of work focuses on 
aggregating constraints into one or more surrogate constraints for which one can show that the 
solution of the original problem and that of the surrogate problem are identical. Early work 
includes that of Balas (1965), Glover (1968), Geoffrion (1969), and Glover (1977); see also the 
survey paper by Rogers et al. (1991). 

We note that the theory of surrogate constraints focuses mainly on integer programs with non¬ 
negative integer coefficients. The seminal work by Mathews (1896) results in an exponential grow 
of the coefficients in the surrogate constraints. Many refinements have been developed, including 
simultaneous and sequential aggregation schemes, to reduce the magnitude of coefficients (Rogers 
et ah, 1991). For mixed-integer programs with real-valued coefficients, however, we are not aware 
of aggregation schemes that guarantee equivalence between the original problem and the surrogate 
problem. In our two-level approach, we employ constraint aggregation as a relaxation technique 
to identify a smaller set of constraints that are active at the optimal solution. This is achieved by 
solving a sequence of MILPs and adding the violated constraints until all constraints in the original 
MILP are satisfied. Furthermore, we take advantage of the LP warm-start to reduce the number 
of MILP re-solves and the computational time of each MILP. 

A large body of literature appeared in the 1990s on bilevel or multilevel mixed integer programs; 
see (Moore and Bard, 1990; Wen and Yang, 1990; Edmunds and Bard, 1992; Hansen et ah, 1992; 
Wen and Huang, 1996) and survey papers (Wen and Hsu, 1991; Vicente and Calamai, 1994). One of 
the main theoretical emphases is on decentralized decision-making from a game-theoretic point of 
view (Wen and Hsu, 1991). Optimality conditions for convex bilevel programs have been established 
in (Vicente and Calamai, 1994). A heuristic-based branch-and-bound method is developed by Wen 
and Yang (1990) and Hansen et al. (1992), and tabu search method is introduced by Wen and Huang 
(1996). Recent years have seen specific application-driven algorithms ranging from infrastructure 
protection planning (Scaparra and Church, 2008) to vulnerability analysis of power grids (Pinar 
et ah, 2010). 

In this context, the multilevel method is specifically referred to the problem formulation. In 
contrast, our two-level method concentrates on algorithmic development for mixed-integer linear 
programs. In particular, our two-level approach builds optimization models of different resolu¬ 
tions from a fine-level problem to a coarse-level problem. We develop a systematic procedure 
that coarsens binary and continuous variables in addition to the aggregation of constraints. Our 
approach allows us to solve large MILPs with more than one million binary and continuous vari¬ 
ables and coupling constraints. Extensive numerical experiments have been conducted to verify the 
efficiency of the developed algorithm. 

Our two-level approach is reminiscent of multi-grid methods in linear algebra and solution of 
partial differential equations. Related work includes that of Gelman and Mandel (1990), who 
discussed a multilevel iterative method for generic optimization problem; Gratton et al. (2008), 
who developed a recursive trust-region method for nonlinear unconstrained problems; and Wen 
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and Goldfarb (2009), who proposed a line search multi-grid approach for nonlinear programs. 

1.2 Our Contributions 

Our contributions are fourfold. First, we develop a systematic two-level approach for two-stage 
MILPs with a large number of binary and continuous variables in the second stage. The coarsening 
of binary variables is done by introducing profiles (vectors) with binary elements. By selecting one 
profile from a small number of candidate profiles, we reduce the number of binary variables by 
orders of magnitude. In addition, we reduce the number of continuous variables by using a convex 
combination of prohles with real elements. We show that the coarsening in variables results in a 
tightening of the original MILP and therefore provides an upper bound on the optimal solution. 

Second, we devise a simple scheme for constraint aggregation that does not require any sparsity 
assumptions on the coefficient matrices. By partitioning constraints into groups and summing up 
constraints in each group, we obtain a relaxation of the original MILP with many fewer constraints. 
We solve a sequence of MILPs by adding any violated constraints until all constraints are satisfied. 
We show that the LP-relaxation of MILPs can be used to warm-start the MILP re-solves and 
significantly reduce the computational effort. 

Third, we implement our two-level approach in the context of the design of cogeneration for 
buildings. The resulting complex MILP model has a large number of coupling constraints over a 
long time horizon. We employ a moving-horizon method to generate valid profiles. We show by 
construction that the generated profiles satisfy coupling constraints in time, hence further reducing 
the number of constraints in the coarsened models. 

Fourth, we verify the efficiency of our algorithm on a variety of examples generated from simu¬ 
lation programs for commercial buildings. While our two-level approach requires solving a sequence 
of coarse MILPs, numerical results indicate that the hrst iterate provides good approximate solu¬ 
tions of the semi-coarse models. Through extensive numerical experiments, we demonstrate the 
scalability of the two-level approach on a rich set of large problems that are beyond the capacity 
of state-of-the-art commercial solvers. 


Outline. The remainder of this paper is organized as follows. In Section 2, we describe the two- 
level approach for the two-stage MILPs. We show that our approach results in a tightening of the 
original problem and provides an upper bound on the optimal solution. In Section 3, we apply our 
approach to a complex MILP model from the cogeneration problem in buildings. We discuss how 
the profiles are generated and selected such that the coarsened models provide feasible solutions 
to (1.1). In Section 4, we demonstrate the effectiveness of our approach on a diverse set of test 
problems. We show that the two-level approach allows us to find good approximate solutions with 
a fraction of computational time compared with solving the full model. In Section 5, we summarize 
our contribution and discuss future extensions. 


2 Two-Level Approach to MILP 

In this section, we derive two models that are formulated in terms of the hrst-stage variables and 
n N second-stage variables. We start by presenting our ideas for variable coarsening, resulting 
in a semi-coarse model that has 0{n) second-stage variables and 0{N) second-stage constraints. 
Next we show how to further coarsen this model by aggregating constraints, resulting in a coarse 
model with both 0{n) second-stage variables and constraints. 
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Figure 1: The first row shows the hourly on/off variables Xt € {0,1} for 48 hours. The second row 
shows the profile representation with two daily profiles Xk £ {0,1}^^. 


2.1 Primal Coarsening 


To coarsen the variables in the model, we introduce a coarsening factor 6 G Z_|_ and group X 
second-stage variables into n = N/5 equal-sized groups of size 6 (the assumption that the groups 
are equal-sized is not critical): 


Xi = 


^ 1)+1 


\ 


for i = 1,..., n. 


( 2 . 2 ) 


We define groups for v and w analogously. These groups correspond to a partition of the second- 
stage variables x,u, and w. 


(2.3) 



/ xl \ 

/ Vl \ 

f Wi \ 

X = 

; , v = 

: L w = 

. 


\ Xn J 

\ Vn / 

\ Wn / 


Next, we introduce (5-profiles, which are fixed parameter values for a group of variables: 

Xk G {0,1}'^ for k = l,...,K. 


(2.4) 


We show in Section 3 how these profiles are generated and selected. Figure 1 illustrates that iV = 48 
hourly on/off variables can be represented by n = 2 daily profiles with length (5 = 24. Now for 
every Xk, we collect a set of Ik operational (5-profiles for V and W: 

Vjk G with LXk < Vjk < UXk for j = 1,..., 4, k = l,...,K, (2.5) 

and 

0 < Wj G for j = 1,..., J, (2.6) 

respectively. 

Given these sets of (5-profiles, we perform a change of variables that replaces the second-stage 
variables (x, v, w) by a reduced set of variables (x, v, w) for the 5-profiles, resulting in a coarsening 
of the second-stage variables. In particular, we set 

K K 

Xi = '^XikXk, '^Xik<l, 4fce{0,1} for i = 1,... ,n. A: = 1,..., W (2.7) 

k=l k=l 






















































6 


Fu Lin, Sven Leyffer, and Todd Munson 


Thus, we have that x 
write V and w as 

G {0,1}'^"', and our 

goal is to create models where Kn <C N. Similarly, we 

K h 

K 4 

Ik 

Vi = 'Yj Yl ^ij^yjk ) 

Y. Y. ^ijk — 1 ! 

Y 0 < Vijk < 1, for z = 1,..., n, (2.8) 

k=l j=l 

k=l j=l 

i=i 

and 


J 

J 


Wi = YjWij^j ’ 

Y^ Wij ^ 1 ) lor i = 1,..., n. (2-9) 


i=i 

j=i 


We note that we do not enforce Wij G {0,1} or Vijk G {0,1}. The reason is that we wish to allow 
more freedom for choosing operational profiles in the second stage as long as they remain feasible. 
This choice also simplifies the coarsened second-stage problem and provides a valid approach if we 
assume that convex combinations of operational (5-profiles Vjk and Wj are feasible, which is often 
the case in practice. 

The partition of variables (x, v, w) implies a partition of the problem matrices B, C, and D of 
(1.1) as 

B = C=[Ci-.----.Cn], D = [Di-.---:Dn], 

where Bi, Ci, and Di G for i = 1,..., n. Next, we define (5-profile matrices as 

X= [Xi : ••• :Xx], V=[Vii:---:VuM-----yKi:---:VKi^], W = [W^ : ■ ■ ■ : Wj], (2.10) 

where Xk, Vjk, and Wj are vectors of length S. We define aggregated (coarse) matrices B, C, and 
D G 

B = [BiX : ■ ■ ■ : BnX], C = [CiV ■.■■■: CnV], D = [D,W : ■ ■ ■ : D^W], (2.11) 

and aggregated cost vectors b,c,d £ as 

b = [b^X :■■■: blxf, c=[clV ■.■■■■. c^Vf, d = [d^W : ■ ■ ■ : d^Wf. ( 2 . 12 ) 

We obtain the following aggregated MILP, which we refer to as the semi-coarse MILP, because 
it has been coarsened in primal variables only: 

minimi^ze a^y iFx -)- (Fv -|- Fw 

y,x,v,w 

subject to Ay Bx -|- Cv + Dw < f 

K 

Xik G {0,1} 

4 (2.13) 

1 ) ^ ^ '^ijk — Xil- , 0 ^ ^ 1 

i=i 

0 < Wij < 1. 

j=l 

The MILP (2.13) has potentially many fewer variables than (1.1), but contains as many constraints 
as the original problem. Before we show how to aggregate constraints in the next section, we finish 
this section by summarizing some of the properties of (2.13). 


Xik < i, 

k=l 
K 4 

< 

k=l j=l 
J 

Wij < 1, 
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Proposition 2.1. Let (x,v,w) be a feasible point of the semi-coarse model (2.13), then it follows 
that the corresponding fine-scale variables 



( K _ \ 

^ ^ ^\k^k 
k=l 


f ^ \ 

/ . xijkVjk 
k,j=l 


- \ 

i=i 

X = 

K 

/ ^ ^nk^k 

Li J 

, V = 

K,Ik 

^njk^jk 

\ k,j=l / 

, W = 

j 

'^WnjWj 
\ i=i / 


are feasible in the original MILP (1.1). 


(2.14) 


Proof. The binary constraint x G {0, l}'^ follows from the representation (2.7) and the definition 
of binary vectors X^. Similarly, the non-negativity of tc is a direct consequence of non-negative 
profiles Wi and non-negative coefficients wa. To show that v is feasible, we start with the definition 
ofVjk in (2.5): 

L-Xk<Vjk<U- Xk. 


Since Vijk > 0, it follows that 


K h K Ik K Ik 

SEE ^ijk^jk — U EE '^ijk^k' 

k=l j=l k=l j=l k=l j=l 


Using Xik = ^ijk^ we have 


K K 

L ^ XikXk <Vi<U'^ XikXk, for i = 1,..., n. 

k=l k=l 


Using the definition of Xi in (2.14), it follows that Lxi < Vi < Uxi, for i = l,...,n, and thus 
Lx < V < Ux. The proof is complete by noting that 


f > Ay + Bx + Cv + Dw = Ay Bx + Cv + Dw, 


where the equality follows from the definition of aggregated matrices (2.11) and the representation 
of the fine-scale variables (2.14). □ 


Next, we show that the semi-coarse model provides an upper bound. 

Proposition 2.2. The semi-coarse model (2.13) is a tightening of the original MILP (1.1), and its 
solution provides an upper bound on (1.1). The two problems are equivalent if the optimal profiles 
from the solution of (1.1) are included in (2.13). 

Proof. Let and optimal value of the original MILP (1.1) and the semi-coarse 

model (2.13), respectively. Since any feasible point of the semi-coarse model (2.13) is a feasible 
point of the original MILP (1.1), it follows that (2.13) is a tightening of (1.1) and thus 


^MILP — ^semi ■ 
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Let {x*, V*, w*) be the optimal solution of MILP (1.1). We now extract optimal profiles {X*, V*, W*) 
corresponding to {x*,v*,w*). Then it follows that there exists a solution of (2.13) such that 

K K Ik J 

4 = < = ee '^ijk^jk^ ^ ^ ^ 1, . . . , U. 

k=l k=l j=l j=l 

Moreover, solution (x*,v*,w*) is feasible in (2.13), and the objective value is the same as solu¬ 
tion [x*,v*,w*). Because (2.13) is a tightening of (1.1), there cannot be a better solution than 
{x*,v*,w*) that we constructed. Hence, {x*,v*,w*) is the optimal solution, and both problems 
(1.1) and (2.13) are equivalent. □ 


2.2 Dual Coarsening 

We next coarsen the constraints by partitioning them into n groups of size 5 and summing over each 
group. We define the rows of A as Af, and we define a set of aggregated constraints by summing 
over the rows of A within each group: 


A : = 


i=l 

S 

EA 

i=l 


T 

■&+i 


<5 

El ^5{n-l)+i 




We aggregate the constraint matrices B, C, and H in a similar way 


B := 


EE^ 
2 = 1 
(5 

E«: 

2=1 


T 

( 5+2 


vEa 

\ i=l 


T 

5(n—l)+i 


C := 


E< 5 ? 

2 = 1 
(5 

E< 5 : 

2 = 1 


-T 

5+2 




\ / ^^S(n-l)+i 

\ i=l 


D := 


Es? 

2=1 

5 

Es: 

2=1 


,T 

5+2 


vE« 

\ i=l 


T 

S{n—l)+i 


Note that we use “bar” notation to denote coarsening in variables and “hat” notation to denote 
coarsening in constraints. We coarsen the right-hand side / analogously and obtain the coarse 
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MILP: 

minimize a?"y + iF x + (Fv + Fw 

v,w,x^y 

subject to Ay + Bx + Cv + Dw < f 

K 

'^Xik<l , Xik G {0,1} 

^=1 

K 4 4 

^ ^ ^ ^ Vijk ^ 1 ) ^ ^ '^ijk — Xik ) 0 ^ Xijk ^ 1 

k=l j=l j=l 

J 

Wij < 1 , 0 < Wij < 1 . 

i=i 


(2.15) 


It follows easily that (2.15) is a relaxation of (2.13), because we have simply aggregated the 
constraints. We can use this fact to develop a simple algorithm that solves (2.13) by solving a 
sequence of tighter relaxations. The main idea is that after solving a relaxation (2.15) we can 
check whether all constraints in (2.13) are satished and add any violated constraints to (2.15). It 
follows easily that this algorithm is finite, because after finitely many constraints have been added 
to (2.15), it is equivalent to (2.13). 

In practice, however, solving a sequence of MILPs (2.15) may not be efficient, because MILPs 
do not warm-start. Instead, we solve the LP-relaxation of the coarse model (2.15) and add vi¬ 
olated constraints until all constraints in the LP-relaxation of the semi-coarse model (2.13) are 
satished. We use the identihed constraints as the initial set of constraints for the MILP coarse 
model iterations. We summarize this procedure in Algorithm I. 


LP Warm-start Phase 


Set I •(— 0; 

Solve the LP-relaxation of the coarse MILP (2.15) to get solution (y^p, x^p, t()°); 
while (ypp, Xpp, ti;^) is not feasible for the LP-relaxation of the semi-coarse MILP (2.13) 

do 


Find constraints that are violated in the LP-relaxation of (2.13) and add them to (2.15); 
Solve the LP-relaxation of (2.15) with the added constraints to get 


Set I — I \ ; 


MILP Phase 


Set A; •(— 0; 

Solve the coarse MILP (2.15) with additional constraints added in the LP Warm-start 

Phase to get solution ( 7 /°, rh'^); 

while (y^,x^,v^,w^) is not feasible for (2.13) do 

Find constraints that are violated in (2.13) and add them to (2.15); 

Solve the coarse MILP (2.15) with the added constraints to get {y^^^ 

Set k = k + 1 ; 


Algorithm 1; Solve the semi-coarse MILP (2.13) via a sequence of coarse MILPs (2.15) using 
the LP-relaxation as warm-start. 
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Proposition 2.3. When Algorithm 1 terminates, the solution , w^) of the coarse model 

(2.15) with added constraints is the solution of the semi-coarse model (2.13). 

Proof. Since {y^, , w^) minimizes the objective function of the semi-coarse model (2.13) and 

satisfies all constraints (2.13), it is the optimal solution of (2.13). □ 


We show in Section 4 that our approach in Algorithm 1 is advantageous; in particular, it 
significantly reduces the number of MILP re-solves, as opposed to the case without LP-relaxation 
as warm-start. We conclude this section by providing a matrix representation of the two-level 
approach using Kronecker products. Recall that for two generic matrices E G p g 

their Kronecker product is defined as 


E®F = 


e\\F ■ ■ ■ e\rnF 

&nlF • • • CjimF 




Using (2.10), we can write 

B = B{I®X), C = C{I®V), D = D{I®W) 

where I is the identity matrix. Similarly, if we let 1 = (1,..., 1) G be the row vector of all ones 
of length <5, then we can express 


B = {I®1)B = {I®1)B{I®X), 
and 


C = {I®1)C = {I®1)C{I®V), b = {I®1)D = {I®1)D{I®W), 

A = {i®i)a, f = {i®i)f. 


3 Application to Cogeneration for Buildings 

In this section, we apply our two-level approach to the cogeneration problem for buildings. Our 
MILP model (A.21) is adapted from models for cogeneration in commercial buildings (Ren and 
Gao, 2010; Pruitt et ah, 2013). In particular, we take linearized models for fuel cells and water 
tank storage from the work of Pruitt et al. (2013), and we penalize on/off operations using switching 
cost as done by Ren and Gao (2010). While the two-level framework described in Section 2 applies 
to generic MILPs (1.1), the MILP model (A.21) for cogeneration entails several complex constraints 
as discussed in Section 3.2. As a result, additional work is required to construct appropriate semi- 
coarse and coarse models. 

We note that our MILP model (A. 21) for the cogeneration problem has the following features 
that are not included in existing models. First, our model has a large number of binary variables, 
on the order of 0(10®), in the second stage. This is orders of magnitude larger than the number of 
binary variables of Siddiqui et al. (2005), Stabler et al. (2009), Ren and Gao (2010), and Pruitt et al. 
(2013). This modeling feature allows considerably more degrees of freedom for on/off operation 
during the life time of new technologies. Second, our model contains three sets of coupling con¬ 
straints that (i) couple first- and second-stage binary variables, (ii) couple second-stage variables for 
different technologies, and (iii) couple second-stage variables over a long time horizon (e.g., 10-20 
years). The number of coupling constraints is on the same order of variables, namely, 0(10®). This 
makes the problem significantly harder because it does not lend itself to decomposition techniques 
such as Lagrangian relaxation. 

In what follows, we describe the main characteristics of the MILP model (Section 3.1), derive the 
semi-coarse and coarse models (Sections 3.2 and 3.3), and discuss profile generation and selection 
(Section 3.4). 
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3.1 MILP Model 

The cogeneration problem consists of two components: the investment decision and the operation 
planning. The investment decision concerns what new technologies to purchase, while the operation 
planning concerns how to dispatch units over a long-term period (e.g., 10-20 years). We formulate 
a two-stage MILP, where the first-stage variables model investment decisions and the second-stage 
variables model equipment operations. Note that both the first and second stages contain integer 
variables. In particular, on/off operations for new technologies in second-stage are made on an 
hourly basis. 

Following Siddiqui et al. (2005); Stabler et al. (2009); Ren and Gao (2010); Pruitt et al. (2013), 
we consider five technologies: batteries, boilers, solid-oxide fuel cells (SOFCs), combined heat and 
power (CHP) SOFCs, and water tank storage. The selection criterion is based on the installation, 
operation, maintenance, fuel consumption, and carbon emission costs, while meeting electricity and 
heating demands. The detailed MILP model is given in (A.21). 

3.2 Semi-Coarse Model 

We begin by introducing the profiles for variables. We denote profiles by upper-case letters with the 
bar notation on top: Xjk G Vo are the on/off profiles for technology j, Pjki G Vp are the production 
profiles associated with Xjk, Uk G Vu are the profiles for power purchased from the utility, Qk G Vq 
are the heat generation profiles from boilers, Bk, BjP G Vb are the power storage and input/output 
profiles for batteries, and Sk,S‘^'^^ G Vs are the heat storage and output profiles for water tanks, 
where Vo,Vp,Vu,Vq,Vb, and Vs are the corresponding set of profiles. We denote the hth. element 
of a profile by Xjk{h) for h G {1,... ,<5}. Time-dependent parameters and of and discount 
factors Yt are concatenated into profiles and G M*^. 

As described in Section 2, the variables in the original MILP model are represented by using 
profiles in the coarsened models. In particular, binary variables Xijt are coarsened to profiles Xjk, 
and continuous variables pijt,ut,bt, st, s°'^^ , and qt are coarsened to profiles Pjki,Uk, Bk, Sk, and 
Sk'^^,Qk > 0, respectively. The coefficients for profiles are denoted by lower-case letters with a 
bar on top; for example, Xijdk G {0, 1} indicates whether profile k is selected on day d for unit 
i of technology j. As modeled in (C.22c), no more than one on/off profile can be selected for 
fixed i,j, and d. This, in conjunction with binary elements of Xjk, results in binary-valued Xijt- 
Similarly, non-negative profiles Pjki, Uk, Bk, Sk, and Qk > 0 and their non-negative coefficients 
in (C.22f), (C.22p), and (C.22q) result in non-negative rtt, st, s°“*, and qt- 

We next turn to constraints that do not couple variables in time, namely, (A.21b), (A.21d), 
(A.21e), (A.21g), (A.21j), (A.21k), (A.21n), and (A.21o). The profile representation for these 
constraints is straightforward; one simply substitutes the summation of profiles as shown in (C.22b), 
(C.22d), (C.22e), (C.22g), (C.22j), (C.22k), (C.22n), and (C.22o), respectively. Note that all 
inequalities in (C.22) are elementwise inequalities. We do not include the production constraint 
(A.21c) in the semi-coarse model. Instead, we require that the production profile Pjki associated 
with the on/off profile Xjk satisfy 

Rf^Xjk < Pjki < nr^Xjk. (3.16) 

This approach is justified because (A.21c) follows by construction due to the convex combination 
of coefficients for production profiles in (C.22f). By an analogous argument, < st in (A.21n) is 
a direct consequence of the following condition imposed on profiles: 

sr^ < Sk. 


( 3 . 17 ) 
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We next discuss how to coarsen variables that are coupled over time periods. The boundary 
conditions (A.21i) and (A.21n) can be expressed as (C.22i) and (C.22m). The maximum power 
purchased constraint (A.21g) can be rewritten as (C.22g), where 1 denotes the vector of all ones 
with length 6. We next turn to switching and storage constraints (A.21f), (A.21h), and (A.211), 
whose prohle representation needs additional notation. Given an on/off profile Xjk G {0,1}*^, we 
can construct the associated switching profile 


Wjkih) = \Xjkih + l)-Xjkih)\. (3.18) 

Now, the switching cost can be included in the objective function in (C.22). To deal with the 
battery storage constraint (A.21h) that couples bt and over the horizon T, we require that the 
pair of profiles {B^, Bjf’) satisfy 

Bk(h + l) = (l-L^)Bk(h) + BlO(h), h = l,...,S-l. (3.19) 

We assign the same coefficient bdk to both sets of profiles Bk and B]^ throughout (C.22). It follows 
that (A.21h) is satisfied except for hours between profiles, namely, t = dd for d < \T)\. Therefore, 
we introduce constraint (C.22h) to guarantee that (A.21h) holds for t = d6 for d < \T)\. We will 
show in Section 3.4 how to generate prohles that satisfy (3.16), (3.17), and (3.19). The heat storage 
constraint (A.211) for j = chp can be expressed as 

^ Sdk { [Ss - (1 - L^)h]Sk + sr^} + Hd+LkEsSk < (Ef/Ef) J2PijdkiPjki, d < \V\, 

kGVs kGVs i,k,l 

where Is G is the identity matrix, Esi G is a matrix with 1 in the {5, 1) entry, and 

Ss G is a Toeplitz matrix with only nonzero elements being 1 at the first upper-subdiagonal. 

3.3 Coarse Model 

We next turn to the coarse model. Note that constraints in the semi-coarse model (C.22) are 
elementwise equalities or inequalities involving daily profiles. We aggregate hourly constraints by 
summing the elements of daily profiles. Using the hat notation on the top to denote the element¬ 
wise summation of profiles, for example, Pjki = Pjki{h): we coarsen 6 hourly constraints into 

a single constraint. For example, the hourly demand constraint (C.22b) is replaced by the daily 
demand constraint (D.23b). Similarly, the maximum hourly power demand (C.22g) is aggregated 
into the maximum daily demand (D.23g). The symmetric breaking constraint (D.23e) has the 
interpretation that the number of times unit i -|- 1 is turned on is no greater than the number 
of times unit i is turned on in day d. The capacity constraints for technologies (D.23d), (D.23j), 
(D.23n), and (D.23o) imply that the daily power/heat output is bounded by the daily capacity of 
the technology units purchased at the hrst-stage. Since (C.22h), (C.22i), and (C.22m) are scalar 
constraints themselves, no aggregation is applied to (D.23h), (D.23i), (D.23m) in the coarse model. 
Also, since the coefficients of profiles are the same for both semi-coarse and coarse models, (D.23c), 
(D.23f), and (D.23p) stay the same as (C.22c), (C.22f), and (C.22q) in the semi-coarse model. 

3.4 Profile Generation and Selection 

Recall that the on/off profiles Xjk must have binary elements and the profiles Pjki, Qk-, Uk^ Bk, 
Sk, and must have non-negative elements. In addition, we impose constraints (3.16), (3.17), 
and (3.19) in formulating the semi-coarse model in Section 3.2. In this section, we discuss how to 
generate valid profiles, and we provide suggestions for profile selections. 
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One approach to generating profiles that satisfy the above constraints is as follows. We solve a 
number of small instances of the original MILP (A.21), take snapshots of the second-stage solutions, 
and extract profiles from these snapshots. For example, consider a four-day MILP (A.21); that is, 
the time horizon in the second-stage problem has only four days. Solving such a four-day model, 
we have four snapshots of daily operation and production; in particular, we have four sets of on/off 
profiles for xtjt G {0,1} and four sets of profiles pijt, ut, h, st, and qjt > 0. We will show in 
Proposition 3.1 that these are valid profiles for the semi-coarse model (C.22). 

The remaining question is how to choose the short-horizon MILPs. Our objective is to generate 
a rich set of profiles that are representative of the optimal solutions for a long-horizon MILP (A.21). 
To this end, we borrow the moving-horizon approach from model predictive control (Garcia et ah, 
1989; Allgower and Zheng, 2000). Given a long-horizon MILP, the idea is to solve MILPs over 
a short window, roll the window forward, and re-solve the new MILP until the window reaches 
the end of the horizon; see Figure 2 for an illustration. We summarize this approach and provide 
additional details in Algorithm 2. 


Data: The parameters for MILP (A.21) with a horizon of D days, and a window of w days 
with w D. 

Result: Profiles {Pjkh Qk, Uk, Bk, Sk, 5^“*} G G and {Xjk,Wjk} G {0,1}^. 

Set A: ^ 0, W = {1,..., (5u;}, TZ = {1,... ,S}; 
while k + w < D do 

Solve MILP (A.21) in the current window t G W; 

Take snapshots of solutions {xip, Pip, Qt, ut, bt, b]^, st, in t G 71; 

Extract profiles {Xjk, Pjki, Qk, Uk, Bk, B]^, Sk, and Wjk G {0,1}'^ using (3.18); 

Roll the window forward by setting A; •(— A; -|- 1; 

if k w < D then Set W = {Sk -|- 1,... ,S{k tc)}, TZ = {Sk -|- 1,..., S{k 1)}; 

if k w = D then Set W = TZ = {Sk -|- 1,..., SH}. 

Algorithm 2: Moving horizon method for profile generation. 


Proposition 3.1. The moving-horizon method described in Algorithm 2 generates non-negative 
profiles {Pjki, Qk, Uk, Bk, Sk, G 1K+ and binary profiles {Xjk, Wjk} G {0, l}”^. Moreover, the 

profile pairs {Xjk, Pjki], {Sk'^^,Sk}, and{Bk,B}P} satisfy (3.16), (3.17), and (3.19), respectively. 

Proof. Non-negativity of the profiles {Pjki, Qk, Uk, Bk, Sk, follows from the fact that the 
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second-stage solutions {pijt, qt, ut, bt, st, of the original MILP model (A.21) are non-negative. 
Similarly, Xjk G {0,1}*^ follows from Xijt G {0,1}; and Wjk, constructed from (3.18), is elementwise 
binary. Since the production variable pijt and the on/off variable xtjt satisfy (A.21c), it follows 
that {Xjk, Pjki} satisfies (3.16). Since the power storage bt and power input/output b]^ satisfy 
(A.21h) and since the heat storage st and heat output satisfy (A.211), we conclude that (3.17) 
and (3.19) follow by construction. □ 

Algorithm 2 potentially generates a huge number of profiles; hence, solving (C.22) and (D.23) 
with all generated profiles may be prohibitive. Instead, we select on/off profiles Xjk that appear 
most frequently in the generated profile pool. Since the aim of the two-level approach is to reduce 
the problem size, it is desired that the number of on/off profiles k is much smaller than the length 
of profiles 6 (e.g., k « <5/10). For the production profiles Pjki associated with each Xjk, we 
choose production profiles that have the minimum or maximum total production Ylh=i 
These extreme profiles provide an envelope of other profiles; thus, their convex combination (C.22f) 
provides a good range of profiles for selection. Similarly, profiles with the minimum and maximum 
sum of absolute values are chosen for the battery storage Bk, the heat storage Sk, the battery 
input/output , and the heat output S'™*. Further, profiles for the power purchased from the 
utility Uk and the heat output from the boiler Qk are uniformly sampled over the period of the 
entire horizon. 

Alternatively, we also employ a k-means clustering algorithm in order to cluster profiles. Given 
a prespecified k number of clusters, this algorithm assigns profiles to one of k clusters defined by 
the centroids (Jain, 2010). Since demand and pricing data for buildings tend to differ significantly 
in winters and summers, we set k = 2, and choose one profile that has the minimum least-squares 
distance to the centroids. In our numerical experiments in Section 4, we apply a /c-means algorithm 
to the boiler output Qk, power purchased from the utility Uk, heat storage output 5™*, and battery 
power output B^ profiles. The clustering approach achieves better performance in the objective 
function than does simple sampling heuristics for profile selection. 

4 Numerical Results and Case Studies 

In this section, we illustrate the effectiveness of our two-level approach using five different building 
examples. We demonstrate that both the semi-coarse and the coarse models allow us to find 
good approximate solutions in a fraction of the time compared with solving the full MILP model. 
The two-level approach also scales to large problems that are beyond the scope of state-of-the-art 
commercial MILP solvers. 

4.1 Generation of Problem Instances 

We use the simulation program Energy Plus 8.2. (Crawley et ah, 2001, 2000) to generate hourly 
electricity and heating demands for typical summer and winter days for five types of buildings 
located in Chicago, Illinois. We scale the daily demands, depending on building types, to generate 
weekly summer and winter demands. The scaling factors for building types are shown in Table 1. 
We generate yearly demand prohles by interpolating between winter and summer weeks. The 
demand of the fth week is given by Di = WiDv/t + (1 — Wi)Dsm, where Uwt and Usm are the weekly 
demand in winter and summer, respectively, and wt is a piecewise linear function shown in Figure 3. 
The resulting yearly demand profiles are perturbed by a zero mean unit variance Gaussian noise 
whose magnitude is 2% proportional to the magnitude of demands. We repeat this process to 
generate demands for multiple years. 
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Building Type 

Mon. 

Tue. 

Wed. 

Thu. 

Fri. 

Sat. 

Sun. 

Office 

1 

1 

1 

1 

1 

1/2 

1/4 

Supermarket 

1/2 

1/2 

1/2 

1/2 

1/2 

1 

1 

Restaurant 

1/4 

1/4 

1/2 

1/2 

1 

1 

1/2 

Hospital 

1 

1 

1 

1 

1 

1 

1 

Retail 

1/4 

1/2 

1/2 

1 

1 

1/4 

1/4 


Table 1: Scaling factors of daily demands in a week for five buildings. For example, the demand of 
the office building on Saturdays is scaled by half of the demand on weekdays. 



Figure 3: Piecewise linear interpolation between winter and summer weeks. The demand of the ith 
week is given by Di = wiD^ + (1 — Wi)Dsm, where D^t is the weekly demand in winter and Dsm 
is the weekly demand in summer. 


Figure 4 shows the electricity demand in summer and winter weeks for five different buildings: 
an office building, a supermarket, a full-service restaurant, a hospital, and a stand-alone retail 
store. In addition to different peak demand patterns, the number of kilowatt-hours varies by orders 
of magnitude for different buildings, from tens of kilowatt-hours for the restaurant to more than 
a thousand kilowatt-hours for the hospital. The diversity of peak demand patterns and the wide 
range of magnitudes show the richness of the example set. As a result, the difficulty of the MILP 
model (A.21) varies significantly for different buildings; see Section 4.3. 

We follow the pricing structure of Stadler et al. (2009). The price for electricity and gas is $0.12 
per kWh and $0,049 per kWh, respectively. The peak demand charge is $14.2 per kW for summer 
months (from June to September) and $11.36 per kW for winter months (from October to May). 
These prices are fixed for each year over the second-stage time horizons. 

4.2 Numerical Experiment Setup 

Numerical experiments are performed on a workstation with 32 GB memory and two Intel E5430 
Xeon 4-core 2.66 GHz GPUs. We implement our algorithms in AMPL (Fourer et ah, 2012) to take 
advantage of AMPL’s compact modeling syntax in profile generation, storage, and management. 
We use GPLEX version 12.6.1.0 as the MIP solver in AMPL. We set a 3-hour time limit and 1% 
relative gap as the stopping criteria for GPLEX. 

4.3 Solutions of the Full Model with Short Second-Stage Horizons 

Table 2 shows the problem size of the full MILP model (A.21) as a function of the number of days 
in the second stage. The number of binary variables, continuous variables, and constraints grows 
linearly with the number of days in the second-stage problem. For the one-year model, there are 
1.05 X 10^ binary variables, 2.71 x 10^ continuous variables, and 6.11 x 10^ constraints. 

Table 3 shows the solutions of small problems for the restaurant example using GPLEX. We 
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(d) Electricity demand for a hospital in winter and summer weeks 



(e) Electricity demand for a stand-alone retail in winter and summer weeks. 
Figure 4: Electricity demand for five buildings in winter and summer weeks. 
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Days 

Binary 

Continuous 

Constraint 

4 

1,152 

2,994 

6,698 

7 

2,016 

5,226 

11,738 

14 

4,032 

10,434 

23,498 

28 

8,064 

20,850 

47,018 

84 

24,192 

62,514 

1.41 • 10® 

364 

1.05 • 10® 

2.71 • 10® 

6.11 • 10® 


Table 2: Problem size of the MILP model (A.21). The number of binary variables, continuous 
variables, and constraints grows linearly with the number of days in the second stage. 


Days 

Time(s) 

Nodes 

LP-iter 

Bat 

Boil 

Chp 

Pow 

Stor 

4 

2 

0 

4.00e-b03 

1 

1 

2 

0 

6 

7 

241 

2807 

4.49e-b05 

2 

1 

0 

1 

0 

14 

1.21e-b03 

10352 

2.73e-b06 

2 

1 

0 

1 

0 

28 

5.69e-b03 

47770 

6.02e-|-06 

2 

1 

0 

1 

0 

84t 

1.06e-b04 

487 

1.31e-b06 

6 

1 

0 

1 

0 


Table 3: Solutions of the MILP model (A.21) for the restaurant example using CPLEX. The time 
in seconds, the number of nodes, and the number of simplex iterations grow exponentially with 
the horizon length in the second stage. The first-stage solutions for batteries (Bat), boilers (Boil), 
CHP-SOFC (CHP), Power SOFC (Pow), and water tank storage (Stor) vary with the problem size. 
^The 84-day model reaches the 3-hour limit, but the relative MIP gap is greater than 18%. 


see the exponential increase of the amount of time, the number of nodes in the branch-and-bound 
trees, and the number of simplex iterations with the number of days. Solving the 84-day example 
is beyond the capabilities of CPLFX on the designated workstation. After reaching the 3-hour 
time limit, the relative MIP gap for the 84-day model is still greater than 18%.^ Our experience 
indicates that it is unlikely that we can solve our MILP over a ten-year time horizon. 

Similar observations can be made for other building examples. In Figure 5, we see that for 
all five buildings, both the number of simplex iterations and the amount of solution time increase 
exponentially with the number of days in the second stage. 

Additional numerical results are summarized in Table 6 in Appendix B. We point out that 
different buildings show different levels of difficulty for the MILP model (A.21). As we see in 
Figure 5, the five buildings differ by orders of magnitude in solution time and the number of simplex 
iterations. Similarly, the number of nodes in the branch-and-bound trees varies significantly over 
different buildings. For example, the 84-day model requires 10,384 nodes for the retail example as 
opposed to 200 nodes for the supermarket example; see Table 6 in Appendix B. 

4.4 Solutions of Semi-Coarse Model with Long Second-Stage Horizons 

In Figure 6, we compare the problem size for the full model (A.21), the semi-coarse model (C.22), 
and the coarse model (D.23) over a 10-year second-stage horizon. We see a roughly 8 times reduction 
in the number of binary variables and 12 times reduction in the number of continuous variables. 
The reason is that for each day with 24 hours, we pick the three most frequent on/off profiles 

^Removing the time limit does not help because then the branch-and-bound tree generated by CPLEX will consume 
all allowable disk space of 100 GB on the workstation. 









Figure 5: Exponential increase of the amount of time and the number of simplex iterations with 
the number of days in the second stage of the full MILP model (A.21). The number of simplex 
iterations drops for the 84-day restaurant example because CPLEX reaches 3-hour time limit. 


for binary variables and two cluster centers from the A;-means clustering algorithm for continuous 
variables. While increasing the number of profiles improves the quality of coarsened models (C.22) 
and (D.23), the resulting computational effort increases significantly. In practice, we find that 
the two-level approach strikes a good balance between solution quality and computational time 
with a small number of profiles (typically 2-3). Note that the large reduction in the number of 
constraints from the full model to the semi-coarse model. This is mainly because we have embedded 
the min/max production constraints (A.21c) and the switching constraints (A.21f) in the profile 
generation; see Section 3.2. In particular, for a 10-year model with 12 fuel cell units, there is a 
reduction of 12 x 87600 x 4 ~ 4.2 x 10® constraints (that is, a 67% of reduction in the number of 
constraints.) 

Table 4 shows the solution information for the semi-coarse model over a 10-year horizon. As 
opposed to the exponential increase for the full model, the solution time, the number of nodes, and 
the number of simplex iterations for the semi-coarse model grow more slowly. Furthermore, the 
first-stage solutions for the semi-coarse model are no longer sensitive to the horizon length. This 
result is in contrast to the variability of first-stage solutions for the full model in Table 3. 

Figure 7 shows that both the number of simplex iterations and the amount of solution time 
scale more slowly for the semi-coarse models. For example, the number of simplex iterations for 
10-year semi-coarse models is on the same order as that of the 84-day full model. The numerical 
results for the semi-coarse model are summarized in Table 7 in Appendix B. 

4.5 Solutions of the Coarse Model with Long Second-Stage Horizons 

As shown in Figure 6, the semi-coarse and coarse models have the same number of binary and 
continuous variables, and the semi-coarse model has about twice as many constraints than does 
the coarse model. Therefore, one expects that both models can be solved with a similar amount of 
computational effort in terms of time and simplex iterations. Since we solve a sequence of coarse 
models, the total computational effort often exceeds that of solving the semi-coarse model. Because 
of the LP warm-start phase in Algorithm 1, however, we find that the first iterate of the coarse 
MILPs provides remarkably good solutions of the semi-coarse model. 
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1 2 4 6 8 10 

Years 


—Full Model Semi-Coarse —Coarse Model 

Figure 6: Problem size for the full MILP model (A.21), the semi-coarse model (C.22), and the coarse 
model (D.23). The number of binary and continuous variables are the same for the semi-coarse 
and coarse models. 


Days 

Time(s) 

Nodes 

LP-iter 

Bat 

Boil 

Chp 

Pow 

Stor 

364 

36.4 

434 

21,700 

0 

1 

1 

1 

1 

728 

157 

673 

57,700 

0 

1 

1 

1 

1 

1,456 

856 

1,540 

1.64-10^ 

0 

1 

1 

1 

1 

2,184 

1,320 

1,737 

2.48 • 10^ 

0 

1 

1 

1 

1 

2,912 

1,520 

1,529 

2.37-10^ 

0 

1 

1 

1 

1 

3,640 

3,350 

2,512 

3.01 • 10^ 

0 

1 

1 

1 

1 


Table 4: Solutions of the semi-coarse model (C.22) for the restaurant example using CPLEX. The 
time in seconds, the number of nodes, and the number of simplex iterations grow much more slowly 
than for the full model in Table 3. The first-stage solutions for batteries (Bat), boilers (Boil), 
CHP-SOFC (CHP), power SOFC (Pow), and water tank storage (Stor) do not change with the 
problem size. 
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Office —SuperMarket —Restaurant Retail —Hospital 


Figure 7: Solution time and nnmber of simplex iterations for the semi-coarse model over a 10-year 
horizon. The growth of computational effort is considerably slower than that of the full model 
shown in Figure 5. 


Iter 

Binary 

Continuous 

Constraint 

Objective Value 

1 

131,040 

254,935 

589,321 

3,797,389.08 

2 

131,040 

254,935 

589,617 

3,797,565.30 

3 

131,040 

254,935 

589,629 

3,797,568.94 


Iter 

Time(s) 

Nodes 

LP-iter 

Bat 

Boil 

Chp 

Pow 

Stor 

I 

7,821.83 

3,620 

1.06 • 10*^ 

0 

1 

1 

1 

1 

2 

1,341.28 

472 

2.38 • 10^ 

0 

1 

1 

1 

1 

3 

1,125.19 

938 

1.47- 10® 

0 

1 

1 

1 

1 


Table 5; Problem size of MILP iterates for a 10-year model of the restaurant example. The hrst 
iterate accounts for 76% of computational time, and the relative gap between the objective valne 
dehned in (4.20) is less than 4.7 x lO”'^. 


To illustrate this point, we show in Table 5 the solntion history of the MILP phase in Algo¬ 
rithm 1. The MILP phase converges in three coarse MILP iterates, and the first iterate requires 
much more computational effort than do other iterates; in particnlar, it acconnts for 76% of total 
computational time. Compared with the solution from the semi-coarse model, the first iterate of 
the coarse MILPs provides a very good solution. The relative gap of the objective value 

= (ObjVal^ejni - ObjValco,,J/ObjVal,emi (4-20) 

between the semi-coarse model and the coarse model at the first iterate is less than 4.7 x 10“^. 

Fignre 8 shows that the relative gap /r for all five building examples over the 10-year horizon is 
no greater than 0.05. This implies that the first iterate of coarse MILPs is a good approximation 
of the semi-coarse model. Numerical results for the first iterate of coarse MILPs are summarized 
in Table 8 in Appendix B. 
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— Office — SuperMarket — Restaurant Retail — Hospital 

Figure 8: The relative gap fj, defined in (4.20) between the first iterate of coarse MILPs and the 
semi-coarse model is no greater than 0.05 for all five buildings over a 10-year horizon. 

4.6 Effect of the LP Warm-Start 

To illustrate the effect of the LP warm-start, we consider the cogeneration problem for the restau¬ 
rant example over a one-year second-stage horizon. The original MILP (A. 21) has 1.05 x 10^ binary 
variables, 2.71 x 10^ continuons variables, and 6.11 x 10^ constraints. Figure 9 shows the num¬ 
ber of constraints and the computation time during the LP warm-start phase. A large number 
of constraints are added in the first few (3-5) iterates, resulting in more computation time at the 
beginning of the algorithm. As the LP warm-start phase progresses, fewer constraints are violated, 
and each LP-relaxation becomes easier to solve because of the warm-start. The amount of time in 
the LP warm-start phase is a fraction of that with one MILP re-solve in the MILP phase. 

We next compare the MILP phase of Algorithm 1 with and without the LP warm-start. Fig¬ 
ure 10 shows the number of MILP re-solves and the compntation time for both cases. Without the 
LP warm-start. Algorithm 1 takes 6 MILP re-solves, as opposed to 2 MILP re-solves with the LP 
warm-start. Moreover, the MILP re-solves with the LP warm-start is up to 20 times faster than 
the MILPs withont the LP warm-start. We note that the improvement of the LP warm-start phase 
is observed in all building examples. 

4.7 First-Stage Solutions 

Since the first-stage solutions are of primary importance, we next compare them from the full 
MILP, the semi-coarse, and the coarse models. In Table 6, we see that the investment decisions 
at the first-stage vary significantly with respect to the horizon length of the second-stage problem. 
For example, the number of water tanks for storage drops from six units in the 4-day model to zero 
units in the 7-day model, and the number of batteries increases from two units in the 28-day model 
to six units in the 84-day model. These numbers imply that the optimal first-stage solutions from 
a short-horizon problem are suboptimal solutions for a long-horizon problem. Therefore, we must 
solve long-horizon MILPs (A.21) as done for the semi-coarse and coarse models. 

In Table 7, the first-stage solutions are much less sensitive to the horizon length at the second 
stage. For example, the first-stage solutions do not change over the 10-year horizon for the super¬ 
market and restaurant examples. For the other three buildings, the £i-norm between first-stage 
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• 10 ^ 




Figure 9: LP warm-start phase for the one-year restaurant example. The number of constraints 
(left panel) and the solution time (right panel) for LP-relaxations of (D.23). 


•lO"^ 




Figure 10: Effect of LP warm-start for the one-year restaurant example. Algorithm 1 with LP 
warm-start phase (o) converges with fewer MILP iterations than without LP warm-start (•). While 
more constraints are added to MILPs with LP warm-start (left panel), each MILP iterate is 5 to 
20 times faster than MILPs without LP warm-start (right panel). 
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solutions of the same building example is no greater than 2 (i.e., at most two-unit difference of 
investment decisions). Similarly, the first-stage solutions at the first iterate of the coarse MILPs 
are insensitive to the horizon length at the second stage. The £i-distance between the first-stage 
solutions of the same building example is no greater than 3. 


5 Conclusions 

We study two-stage MILPs with more than 1 million binary variables in the second-stage problem. 
We develop a two-level approach by constructing semi-coarse models (coarsened with respect to 
variables) and coarse models (coarsened in both variables and constraints). We show that the semi- 
coarse model is guaranteed to provide a feasible solution of the original MILP and hence resnlts in 
an upper bound on the optimal solution. We solve a sequence of coarse MILPs with aggregated 
constraints that converges to the same upper bound with a finite number of steps. Furthermore, 
we take advantage of the LP warm-start to reduce the number of MILP re-solves. 

We apply our approach to the cogeneration problem for commercial buildings. We demonstrate 
the effectiveness of the two-level approach using building examples with simulation data. In par- 
ticnlar, the two-level approach allows us to obtain good approximate solutions at a fraction of the 
time required for solving the original MILPs. Furthermore, we show that the two-level approach 
scales to large problems that are beyond the capacity of state-of-the-art commercial solvers. 

A number of extensions are of interest in our future work. First, how should one add new 
profiles in the coarsened models? Currently, we select a number of prespecified profiles and fix 
them throughout the MILP re-solves. Since the solution qnality of coarsened models is determined 
by profiles, it is of interest to add new profiles dynamically that are most promising in reducing the 
objective valne. It seems that the reduced cost associated with profiles can be obtained by solving 
appropriate pricing problems like those in the column generation approach (Barnhart et ah, 1998). 
Further, it is also of interest to consider dynamic aggregations of constraints as new profiles are 
included. For set-partitioning constraints, a similar approach has been developed in the column 
generation framework (Elhallaoui et ah, 2005). We plan to explore this line of research in future 
work. 

Second, how can one obtain lower bounds on the optimal solution in the two-level framework? 
Since the coarse model is a relaxation of the semi-coarse model, which itself is a tightening of 
the full model, the coarse model is not a relaxation of the original MILP. While relaxing binary 
variables in the original MILP resnlts in a lower bound, our numerical results indicate a sizable gap 
between this lower bound and the optimal solution. It is an open problem how to generate lower 
bounds using appropriate coarse MILPs. 

Furthermore, advanced multilevel approaches for PDFs typically require several sweeps of 
fine/coarse levels in order to achieve accurate solutions (Falgout, 2006). It may be beneficial 
to iterate between solving coarse models to optimality and solving “partially” the fine model until 
a prespecified accuracy for the solution is achieved. We intend to extend our two-level framework 
in this direction. 
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A Full Model for the Cogeneration Problem 

For the sake of completeness, we describe here our cogeneration model. 

Index Sets. We use calligraphic upper-case letters for all sets. The set of technologies: batteries, 
boilers, CHP-SOFC, power SOFC, and water tank storage is denoted by J. We use batt, boil, 
pow, chp, and stor for shortcuts. The subset Jg = {power SOFC, CHP-SOFC} denotes SOFC 
generation technologies that require on/off operations at the second stage. Each technology j ^ Jg 
has a number of identical units, and the index set is denoted by lAj. Finally, T is the index set of 
time, M is the index set of months, and Tm C T is the index set of time for each month m € Ai. 

Variables. We use lower-case letters to denote variables. The first-stage variables in our model 
are the number of units of technology j G J, which we denote by G Z_|_. All other variables model 
operation of the installed units, and they are our second-stage variables. We use the subscript t 
to indicate the time period: bt and b]^ are the power storage and input/output for batteries, 
respectively; s*, s™, and are the heat storage, input, and output for water tanks, respectively; 
Pjt and Qjt > 0 are the power and heat output of technology j, respectively; ut > 0 is the power 
purchased from the utility and > 0 is the maximum power purchased in month m; Xijt G {0,1} 
indicates whether unit i of technology j operates in time period t; and Wijt G [0,1] is an auxiliary 
variable to indicate whether unit i of technology j switches from t to t + 1. Note that we do not 
impose binary constraint on Wijt because it takes a binary value at the solution due to the switching 
constraint (A.21f) and the binary constraint Xijt G {0,1}. We use the convention t G T, j £ J, 
* e ^j, m G Ai unless explicitly mentioned. 

Parameters. We use upper-case letters to denote parameters and constants: Cj is the capital 
and installation cost of technology j; H is the number of hours in the lifetime of technology (e.g., 
H = 87600 for a ten-year model); T is the number of hours in the problem horizon; Y is the 
hourly discount factor with 3% annual interest rate (i.e., Y = 1 — 0.03/8760); Mj is the operation 
and maintenance cost of technology j; Pt and Gt are the price for electricity and natural gas, 
respectively; is the peak demand price for power from the utility in month m; Wj is the cost 

of switching a unit of technology j on or off; 5} is the thermal capacity per unit for technology 

j; is the power capacity per unit for batteries; Df and of are the power and heat demand 
in period t, respectively; 7?™™ and are the minimum and maximum power output when a 

unit of technology j is turned on, respectively; and are the power and thermal efficiency of 
technology j, respectively; and are the average power loss in battery and heat loss in water 
tank storage, respectively. 

Cost Function. We now formulate the complete MILP model for cogeneration problem in (A.21). 
The objective function consists of two parts: the capital and installation cost in the first-stage and 
the operation cost in the second stage. The operation cost includes the peak power usage, operation 
and maintenance, switching units, gas, and purchased power costs. We discount the second-stage 
costs with an hourly discount factor based on a 3% annual interest rate; hence Y = 1 — 0.03/8760. 
Therefore, the cost incurred in hour t is multiplied by the discount factor The peak demand 
charge at the end of month m is discounted with F*’". To compare cost over different time horizons 

^For comparison, the annual discount factor is rs 97.044%, which is consistent with the 3% interest rate; 

and a ten-year discount factor is ~ 74.082%. 
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T, we compute the average hourly cost and multiply the result with the lifetime of technologies, 
H, resulting in the factor H/T in the second-stage cost in (A.21). 


minimize 


subject to 


i&J m&M 

+ E^'[Ew«‘ + E WjWijt) + Gtqt + PtUt] I 

teT jejg i&Uj 


yj G Z, 0 <yj < \l(j\ 


(A.21a) 

Pjt + ut- bl^ > D[ 


(A.21b) 

j&Js 



< Pijt < Pjt = Y PLt 

j G Jg 

(A.21c) 

ieUj 



Xijt G {0,1}, ^ ^ Xijt E yj 

j G Jg 

(A.21d) 

ieUj 



^{i+l)jt — ^ijt 

j G Jg, 

i < \Uj\ 



(A.21e) 

^ij{t+l) ^ijt E Wijt, Xijt ^ij{t+l) — ^ijti 0 E Wijt E 1) 

j G Jg, 

t<\T\ 



(A.21f) 

0<ut< uZ^ 

t G I'm 

(A.21g) 

bt+i = (1 - LP)bt + bf 

t<\T\ 

(A.21h) 

II 


(A.21i) 

0 < 


(A.21j) 

sT" + qt> Df 


(A.21k) 

st+i - (1 - LQ)st + sr < (i^i/^fhp)Pchp,i 

t<\T\ 

(A.211) 

= S\T\ 


(A.21m) 

<St< 5,^torystor 


(A.21n) 

qt E ‘S'boiiyboii 


(A.21o) 

bt,pjt,qt,st,s°'^*‘,ut > 0 


(A.21p) 


Constraints. The constraints in the MILP model (A.21) can be divided into three groups. The 
first group of constraints couples first- and second-stage variables; in particular, the number of 
on-units for fuel cells and the capacity of battery, storage, and boiler are bounded by the num¬ 
ber of units purchased in the first-stage, as described in (A.21d), (A.21j), (A.21n), and (A.21o), 
respectively. The second group of constraints couples variables across technologies; in particular, 
technologies are coordinated to meet the power demand (A.21b) and the heat demand (A.21k). 
We note that dualizing these constraints results in decoupled variables for each technology; see, 
for example, Fisher (1981) and Bard (1988). However, MILP (A.21) does not lend itself to this 
Lagrangian relaxation approach because of coupling constraints in the third group. This group of 
constraints couples variables over time periods in the second stage. In particular, the battery and 
storage constraints (A.21h) and (A.211) link variables over the entire horizon. In addition, con¬ 
straint (A.21g), which models the maximum power purchased from the utility in each months, links 
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variables over time instances in each month. Moreover, constraint (A.21f), which models on/off 
switches of the generation units, couples binary variables Xijt over two immediate time instances. 
Boundary conditions (A.21i) and (A.21m) for battery and storage couple variables at both ends of 
the entire horizon. 

We conclude this subsection by explaining the rest of the constraints in (A.21). The number of 
purchased units is upper bounded in (A.21a). Since all units of the same technology are identical, 
we introduce a symmetry breaking constraint (A.21e) that states that if unit i + 1 is on, then unit i 
must be on as well. When unit i of technology j is turned on, the power generated pijt is bounded 
by the minimum and maximum capacity in (A.21c). The non-negativity constraint for variables is 
given by (A.21p). 

B Tables of Numerical Results 
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Table 6: Problem size and AMPL/CPLEX solution information for the cogeneration problem (A.21) for five buildings: 1-Office, 2-Supermarket, 
3-Restaurant, 4-Retail, 5-Hospital. Model instances that reach 3-hour time limit are indicated by MSG=1. 
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-10^ 

0 

5 

1.55 

• 10-2 

2,912 

1.05 • 103 

2.19-103 

1.48- 

10® 

5,686 

6.29-103 

4 

6 

5 

7 

2 

10,500 

6.59 

- 10^ 

1 

5 

3.64 

• 10-2 

3,640 

1.31 • 103 

2.73 - 103 

1.85 - 

10® 

2,966 

7.87-103 

3 

6 

5 

7 

2 

10,800 

6.55 

- 10^ 

1 


Table 7: Problem size and AMPL/CPLEX solution information for the semi-coarse model (C.22) for five buildings: l-OfRce, 2-Supermarket, 3- 
Restaurant, 4-Retail, 5-Hospital. Model instances that reach 3-hour time limit are indicated by MSG=1. 
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Building 

Year 

Binary 

Continuous 

Constraint 

Nodes 

LP-iter 

Time 

Bat 

Boil 

Chp 

Pow 

Stor 

ObjVal 

MSG 

1 

1 

13,104 

25,498 

67,251 

550 

34,876 

35.09 

0 

3 

5 

7 

2 

2.8 

10^ 

0 

1 

2 

26,208 

50,991 

1.34 - 10® 

2,507 

1.53 - 10® 

324.85 

0 

3 

5 

7 

2 

2.76 

-10^ 

0 

1 

4 

52,416 

1.02 • 10® 

2.73 - 10® 

504 

1.29 - 10® 

415.87 

0 

3 

5 

7 

2 

2.68 

- 10^ 

0 

1 

6 

78,624 

1.53-10® 

4.13-10® 

6,166 

4.52 - 10® 

3,104.26 

0 

3 

5 

7 

2 

2.61 

- 10^ 

0 

1 

8 

1.05- lO*^ 

2.04 • 10® 

5.36-10® 

8,482 

6.13 - 10® 

6,312.24 

0 

3 

5 

7 

2 

2.53 

-10" 

0 

1 

10 

1.31 • 10*^ 

2.55 • 10® 

6.85 - 10® 

7,864 

7.42 - 10® 

10,453.39 

0 

3 

5 

7 

2 

2.48 

-10^ 

1 

2 

1 

13,104 

25,498 

68,345 

1,410 

49,815 

96.47 

0 

3 

5 

4 

2 

1.84 

- 10^ 

0 

2 

2 

26,208 

50,991 

1.4-10® 

3,774 

1.21 - 10® 

464.67 

0 

3 

5 

5 

2 

1.8 

10^ 

0 

2 

4 

52,416 

1.02 • 10® 

2.95 - 10® 

11,443 

3.05 - 10® 

3,019.47 

0 

3 

5 

4 

2 

1.75 

- 10^ 

0 

2 

6 

78,624 

1.53-10® 

4.37-10® 

10,343 

3.63 - 10® 

3,900.92 

0 

3 

5 

4 

2 

1.7- 

10^ 

0 

2 

8 

1.05 • 10*^ 

2.04 - 10® 

5.87-10® 

9,474 

4.39 - 10® 

7,293.15 

0 

3 

5 

4 

2 

1.66 

- 10^ 

0 

2 

10 

1.31 • 10^5 

2.55 - 10® 

7.37 - 10® 

11,718 

5.6 - 10® 

10,683.5 

0 

3 

5 

7 

2 

1.62 

- 10^ 

1 

3 

1 

13,104 

25,498 

55,782 

520 

38,780 

46.85 

0 

1 

1 

1 

1 

4.33 

- 10® 

0 

3 

2 

26,208 

50,991 

1.12 - 10® 

471 

64,663 

106.27 

0 

1 

1 

1 

1 

4.26 

- 10® 

0 

3 

4 

52,416 

1.02 - 10® 

2.32 - 10® 

3,474 

2.2 - 10® 

1,102.5 

0 

1 

1 

1 

1 

4.13 

- 10® 

0 

3 

6 

78,624 

1.53-10® 

3.56-10® 

2,519 

2.06 - 10® 

1,215.44 

0 

1 

1 

1 

1 

4.03 

- 10® 

0 

3 

8 

1.05 • 10*^ 

2.04 - 10® 

4.82 - 10® 

10,529 

6.58 - 10® 

7,214.81 

0 

1 

1 

1 

1 

3.9 

10® 

0 

3 

10 

1.31 • lO'^ 

2.55 - 10® 

5.89 - 10® 

3,620 

1.06 - 10® 

7,821.83 

0 

1 

1 

1 

1 

3.8 

10® 

0 

4 

1 

13,104 

22,950 

56,465 

1,600 

1.66 - 10® 

256.57 

0 

1 

2 

1 

1 

6.14 

- 10® 

0 

4 

2 

26,208 

45,895 

1.15 - 10® 

2,508 

4.11 - 10® 

721.43 

0 

1 

2 

1 

1 

6.05 

- 10® 

0 

4 

4 

52,416 

91,785 

2.29-10® 

38,110 

2.39 - 10® 

7,938.05 

0 

1 

2 

1 

1 

5.88 

- 10® 

0 

4 

6 

78,624 

1.38-10® 

3.37-10® 

45,080 

1.15 - 10® 

10,705.24 

0 

1 

1 

1 

1 

5.71 

- 10® 

1 

4 

8 

1.05- lO'^ 

1.84 - 10® 

4.61 - 10® 

25,554 

6.42 - 10® 

9,382.93 

0 

1 

1 

1 

1 

5.55 

- 10® 

1 

4 

10 

1.31 • 10*^ 

2.29-10® 

5.81 - 10® 

15,413 

6.33 - 10® 

9,669.02 

0 

1 

1 

1 

1 

5.4- 

10® 

1 

5 

1 

13,104 

27,318 

68,979 

1,759 

90,264 

178.67 

4 

7 

5 

7 

2 

7.26 

- 10^ 

0 

5 

2 

26,208 

54,631 

1.39-10® 

3,169 

1.73 - 10® 

761.49 

4 

6 

5 

7 

2 

7.14 

- 10^ 

0 

5 

4 

52,416 

1.09-10® 

2.78-10® 

8,328 

3.76 - 10® 

3,633.01 

3 

6 

5 

7 

2 

6.95 

-10" 

0 

5 

6 

78,624 

1.64 - 10® 

4.16-10® 

10,835 

5.36 - 10® 

6,104.61 

4 

6 

5 

7 

2 

6.75 

-10^ 

0 

5 

8 

1.05 • 10*^ 

2.19-10® 

5.55 - 10® 

10,707 

7.04 - 10® 

10,826.62 

4 

6 

5 

7 

2 

6.58 

- 10^ 

1 

5 

10 

1.31 • 10^5 

2.73 - 10® 

6.91 - 10® 

7,376 

7.96 - 10® 

10,227.72 

4 

6 

5 

7 

2 

6.42 

- 10^ 

1 


Table 8: Problem size and AMPL/CPLEX solution information for the first iteration of the coarse model (C.22) for five buildings: 1-Office, 2- 
Supermarket, 3-Restaurant, 4-Retail, 5-Hospital. Model instances that reach 3-hour time limit are indicated by MSG=1. 
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C Semi-Coarse Model for the Cogeneration Problem 


minimize ^ Cjyj + y{ + Y PijdkiMjPjki 

“1“ ^ ^ ^ ^ QdkQk “1“ ^ ^ '^dk^k'] J" 

dA^Vq dA^Vu 

subject to ifj E Z, 0 < < \Uj\ 3 ^ J 

Y, PijdklPjkl + ^ ^ hdkBf > 

j,i,k,l k&Vu k&Vb 

Xijdk £ { 0 , 1 }, ^ ^ Xijdk ^ 1 

k&Vo 

^ ^ ^ijdkYjk Y Uj^ 
i,k 


Y1 ^{i+Ljdk- 

^jk — ^ ^ ^ijdk^jk 

i < \Uj\ 

k^'Po 

k^Y 0 


Pijdk ^ [0, 1], 

^ ^ Pijdkl — ^ijdki ^ ^ ^ ^ Pijdkl Y 1 



l^Yp k^Yo l^Yp 


0 < ^ UdkUk < 

d G 2?m 

k&Vu 



Yj ^{d+l)kBj 

,{!)= Ybdk[{l-YBk{S) + B]P{S)] 

d< \V\ 

k^Pb 

k&Vb 


Y. bikBk{l) 



kGPb 

k&Vb 


Y. bdkBk < 

B ybattl 


keVb 



Y, sdkST^ + Y ^ ^d 


k^Ys 

k&Vq 


Y. ^dk [5'5 - 

{l-LP)Is]Sk+ Y Hd+iikEsSk 

d< \V\ 

k^Ys 

k&Vs 


+ Y MkSr^ < {EYE!’)YPLdklPjkl 

j = chp 

k^Ys 

i,k3 


Y ^^kSk{l) 

= Y/ ^Ipk^ki^) 


k^Ys 

k^Ys 


Y ~^dkSY < 5’s^tori/storl 


k^Ys 



^ ^ QdkQk Y 

‘S'boil2/boill 


kSYq 



'^dkj 

Qdk £ [0)1] 


^ ^ ^dk — I 5 

^ ^ i^dk — 1) ^ ^ ^ 1, ^ ^ Qdk — !• 


k^Yu 

fcePt fcePs fcePg 



(C.22a) 

(C.22b) 

(C.22c) 

(C.22d) 

(C.22e) 

(C.22f) 

(C.22g) 

(C.22h) 

(C.22i) 

(C.22j) 

(C.22k) 

(C.221) 

(C.22m) 

(C.22n) 

(C.22o) 

(C.22p) 

(C.22q) 
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D Coarse Model for the Cogeneration Problem 


minimize 


subject to 


E <^01 + f { E +EE 


jkl 


H“ ^ ^ ^ ^ ‘^dkUk H“ ^ ^ 

dj^i^k d^k^Vu d,kGVq 

Vj G Z, 0 < < im 

Z PijdklPjkl + '^ UdkUk - hdkBf > 
j,i,k,l k£Pu k^Vb 

Xijdk £ {0) 1}) ^ ^ Pijdk E 1 

k&Vo 

^ ^ ^ijdk-^jk “Ti S ' Vj 


i,k 


^ ^ Xijdkl^jk P ^ ^ 


jfc 


Pijdk £ [ 0 ) 1 ]) ^ ^ Pijdkl ^ijdkj ^ ^ ^ ^ Pijdkl — 1 
l&Pp k&Vo lePp 


0 < ^ UdkUk < S 


U„ 


k£Pu 


E hd+i)kMi) = E + Bf{ 5 )) 

k&Pb kePb 

E ^ikBk{l) = E hv\kBk{S) 

kePb kePb 

^ ^ bdkBk ^ 5 ■ B 2 /batt 
k^'Pb 

E + E ^ 

kePs k&Pq 

E SdkiL^Bk - Skil)) + E '®(d+l)fc'^fc(^) 


k^'P s 


k^Ps 


+ E < {EflEP)Y,PTdkiPjki 


k^P s 


i^k^l 


E ^ikSk{i) = E ^i»ifc'S'fe('5) 

kePs k&Ps 

E«<^fc^r*<<^-5’s"torystor 

k^Ps 

^ ^ QdkQk — ^ ‘ 5*j^Qjj^}2,oil 


k^Pn 


^dki ^dki ^dkt Qdk ^ [^d] 

^ ^ Udk E 1) ^ ^ bdk E 1) ^ ^ Sdk E 1; ^ ^ Qdk E 1- 


kGPu 


k&Pb 


kePs 


k&Pa 


j^J 


i < \Uj 


d G Pm 
d< \V\ 


d< \V\ 
j = chp 


(D.23a) 

(D.23b) 

(D.23c) 

(D.23d) 

(D.23e) 

(D.23f) 

(D.23g) 

(D.23h) 

(D.23i) 

(D.23j) 

(D.23k) 

(D.231) 

(D.23m) 

(D.23n) 

(D.23o) 

(D.23p) 



32 


Fu Lin, Sven Leyffer, and Todd Munson 


References 

Alguacil, N., Motto, A. L., and Conejo, A. J. (2003). Transmission expansion planning: a mixed- 
integer LP approach. IEEE Transactions on Power Systems, 18(3):1070-1077. 

Allgower, F. and Zheng, A. (2000). Nonlinear Model Predictive Control, volume 26. Birkhauser 
Basel. 

Balas, E. (1965). An additive algorithm for solving linear programs with zero-one variables. Oper¬ 
ations Research, 13(4):517-546. 

Bard, J. F. (1988). Short-term scheduling of thermal-electric generators using Lagrangian relax¬ 
ation. Operations Research, 36(5):756-766. 

Barnhart, C., Johnson, E. L., Nemhauser, G. L., Savelsbergh, M. W., and Vance, P. H. (1998). 
Branch-and-price: Column generation for solving huge integer programs. Operations Research, 
46(3):316-329. 

Birge, J. R. (1985). Aggregation bounds in stochastic linear programming. Mathematical Program¬ 
ming, 31(1):25-41. 

Birge, J. R. and Louveaux, F. (2011). Introduction to Stochastic Programming. Springer Science 
& Business Media. 

Clay, R. L. and Grossmann, I. E. (1997). A disaggregation algorithm for the optimization of 
stochastic planning models. Computers & Chemical Engineering, 21(7):751~774. 

Crawley, D. B., Lawrie, L. K., Winkelmann, F. C., et al. (2001). Energyplus: creating a new- 
generation building energy simulation program. Energy and Buildings, 33(4):319-331. 

Crawley, D. B., Pedersen, C. O., Lawrie, L. K., and Winkelmann, F. C. (2000). Energyplus: energy 
simulation program. ASHRAE Journal, 42:49-56. 

Edmunds, T. A. and Bard, J. F. (1992). An algorithm for the mixed-integer nonlinear bilevel 
programming problem. Annals of Operations Research, 34(1):149-162. 

Elhallaoui, I., Villeneuve, D., Soumis, F., and Desaulniers, G. (2005). Dynamic aggregation of 
set-partitioning constraints in column generation. Operations Research, 53(4):632-645. 

Falgout, R. D. (2006). An introduction to algebraic multigrid. Computing in science & engineering, 
8:24. 

Fisher, M. L. (1981). The Lagrangian relaxation method for solving integer programming problems. 
Management science, 27(1):1~18. 

Fourer, R., Gay, D. M., and Kernighan, B. W. (2012). AMPL: A Modeling Language for Mathe¬ 
matical Programming. Cengage Learning. 

Garcia, C. E., Prett, D. M., and Morari, M. (1989). Model predictive control: Theory and practice 
- a survey. Automatica, 25(3):335-348. 

Gelman, E. and Mandel, J. (1990). On multilevel iterative methods for optimization problems. 
Mathematical Programming, 48(1-3): 1-17. 



Two-Level MILP 


33 


Geoffrion, A. M. (1969). An improved implicit enumeration approach for integer programming. 
Operations Research, 17(3):437-454. 

Geoffrion, A. M. (1974). Lagrangean relaxation for integer programming. Springer. 

Glover, F. (1968). Surrogate constraints. Operations Research, 16(4):741-749. 

Glover, F. (1977). Heuristics for integer programming using surrogate constraints. Decision Sci¬ 
ences, 8(1);156-166. 

Gratton, S., Sartenaer, A., and Toint, P. L. (2008). Recursive trust-region methods for multiscale 
nonlinear optimization. SIAM Journal on Optimization, 19(l):414-444. 

Hansen, P., Jaumard, B., and Savard, G. (1992). New branch-and-bound rules for linear bilevel 
programming. SIAM Journal on Scientific and Statistical Computing, 13(5):1194-1217. 

Hedman, K. W., Ferris, M. C., O’Neill, R. P., Fisher, E. B., and Oren, S. S. (2010). Co-optimization 
of generation unit commitment and transmission switching with N-1 reliability. Power Systems, 
IEEE Transactions on, 25(2):1052-1063. 

Jain, A. K. (2010). Data clustering: 50 years beyond K-means. Pattern recognition letters, 
31(8):651-666. 

Mathews, G. B. (1896). On the partition of numbers. Proceedings of the London Mathematical 
Society, l(l):486-490. 

Moore, J. T. and Bard, J. F. (1990). The mixed integer linear bilevel programming problem. 
Operations Research, 38(5):911-921. 

Munoz, F. D., Hobbs, B., and Watson, J.-P. (2014). New bounding and decomposition approaches 
for MILP investment problems: Multi-area transmission and generation planning under policy 
constraints. Technical report, Sandia National Laboratories. 

Munoz, F. D., Sauma, E. E., and Hobbs, B. E. (2013). Approximations in power transmission 
planning: implications for the cost and performance of renewable portfolio standards. Journal 
of Regulatory Economics, 43(3):305-338. 

Pinar, A., Meza, J., Donde, V., and Lesieutre, B. (2010). Optimization strategies for the vulnera¬ 
bility analysis of the electric power grid. SIAM Journal on Optimization, 20(4): 1786-1810. 

Pruitt, K. A., Braun, R. J., and Newman, A. M. (2013). Evaluating shortfalls in mixed-integer 
programming approaches for the optimal design and dispatch of distributed generation systems. 
Applied Energy, 102:386-398. 

Pruitt, K. A., Leyffer, S., Newman, A. M., and Braun, R. J. (2014). A mixed-integer nonlinear 
program for the optimal design and dispatch of distributed generation systems. Optimization 
and Engineering, 15:167-197. 

Ren, H. and Gao, W. (2010). A MILP model for integrated plan and evaluation of distributed 
energy systems. Applied Energy, 87(3): 1001-1014. 

Rogers, D. E., Plante, R. D., Wong, R. T., and Evans, J. R. (1991). Aggregation and disaggregation 
techniques and methodology in optimization. Operations Research, 39(4):553-582. 



34 


Fu Lin, Sven Leyffer, and Todd Munson 


Scaparra, M. P. and Church, R. L. (2008). A bilevel mixed-integer program for critical infrastructure 
protection planning. Computers & Operations Research, 35(6):1905-1923. 

Siddiqui, A. S., Marnay, C., Bailey, O., and Hamachi LaCommare, K. (2005). Optimal selec¬ 
tion of on-site generation with combined heat and power applications. International Journal of 
Distributed Energy Resourees, l(l):33-62. 

Stadler, M., Marnay, C., Siddiqui, A., Lai, J., Coffey, B., and Aki, H. (2009). Effect of heat 
and electricity storage and reliability on microgrid viability: a study of commercial buildings 
in California and New York states. Lawrenee Berkeley National Laboratory Teehnical Report 
LBNL-1334E. 

Vicente, L. N. and Calamai, P. H. (1994). Bilevel and multilevel programming; A bibliography 
review. Journal of Global Optimization, 5(3):291~306. 

Wen, U. P. and Hsu, S. T. (1991). Linear bi-level programming problems-a review. Journal of the 
Operational Researeh Soeiety, pages 125-133. 

Wen, U. P. and Huang, A. D. (1996). A simple tabu search method to solve the mixed-integer 
linear bilevel programming problem. European Journal of Operational Research, 88(3):563-571. 

Wen, U. P. and Yang, Y. H. (1990). Algorithms for solving the mixed integer two-level linear 
programming problem. Computers & Operations Research, 17(2):133-142. 

Wen, Z. and Goldfarb, D. (2009). A line search multigrid method for large-scale nonlinear opti¬ 
mization. SIAM Journal on Optimization, 20(3):1478-1503. 

Zipkin, P. H. (1980a). Bounds for row-aggregation in linear programming. Operations Research, 
28(4):903-916. 

Zipkin, P. H. (1980b). Bounds on the effect of aggregating variables in linear programs. Operations 
Research, 28(2):403-418. 


The submitted manuscript has been created by the UChicago Argonne, LLC, Operator of Argonne National Laboratory (“Ar- 
gonne”) under Contract No. DE-AC02-06CH11357 with the U.S. Department of Energy. The U.S. Government retains for 
itself, and others acting on its behalf, a paid-up, nonexclusive, irrevocable worldwide license in said article to reproduce, 
prepare derivative works, distribute copies to the public, and perform publicly and display publicly, by or on behalf of the 
Government. 




